---
title: "Longletter - Commission Responsiveness"
subtitle: "Tables and Figures"
author: "Christoph Ivanusch, Leonce Röth"
date: "`r Sys.Date()`"
output: 
  pdf_document:
    keep_md: true
    keep_tex: true
link-citations: yes
urlcolor: red
toc: yes
header-includes:
  - \usepackage{float}
  - \usepackage{subfig}
  - \captionsetup{belowskip=0pt,aboveskip=4pt}
  - \usepackage{setspace}
  - \usepackage{lscape}
  - \usepackage{siunitx}
---

```{r, setup, include=FALSE, echo=FALSE, warning=FALSE, message=FALSE}

# chunk option defaults
knitr::opts_chunk$set(echo = FALSE, message = FALSE)

# packages
library(tidyverse)
library(dplyr)
library(zoo)
library(readxl)
library(modelsummary)
library(vars)
library(ggplot2)
library(kableExtra)
library(gridExtra)
library(readxl)

```

\newpage

```{r}
# Set the working directory
setwd("S:/DFG/Backup/Backup - Recent investigations/Talking time/Long_Letter/Data")

# Confirm the working directory is set correctly
getwd()

main_data <- readRDS("main_data.RDS") %>%
  filter(month <= "2023-12-01")

# Read the Excel file
main_data2 <- read_excel("main_data_new.xlsx") %>%
 filter(month <= "2023-12-01")

```

```{r}
library(dplyr)

```

```{r}
# Step 2: Perform the join and update variables
main_data_updated <- dplyr::left_join(
  main_data,
  main_data2 %>% dplyr::select(month, monthly_resp_f_count, monthly_tp_h_count, total_response_count),
  by = "month"
) %>%
  dplyr::mutate(
    monthly_resp_f_count = dplyr::coalesce(monthly_resp_f_count.y, monthly_resp_f_count.x),
    monthly_tp_h_count = dplyr::coalesce(monthly_tp_h_count.y, monthly_tp_h_count.x),
    total_response_count = dplyr::coalesce(total_response_count.y, total_response_count.x)
  ) %>%
  dplyr::select(-monthly_resp_f_count.x, -monthly_tp_h_count.x, -total_response_count.x)

# View the updated dataset
main_data_updated

```

```{r}
# Rename the data frames
main_data_old <- main_data  # Renaming original data to main_data_old
main_data <- main_data_updated  # Renaming the updated data to main_data

```

\newpage

# Transform data

```{r, create time lag variables, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# create time lag variables
main_data <- main_data %>%
  mutate(monthly_resp_f_count_lag1 = dplyr::lag(monthly_resp_f_count, n = 1)) %>%
  mutate(monthly_tp_h_count_lag1 = dplyr::lag(monthly_tp_h_count, n = 1)) %>%
  mutate(asylum_total_lag1 = dplyr::lag(asylum_total, n = 1)) %>%
  mutate(asylum_total_scaled_lag1 = dplyr::lag(asylum_total_scaled, n = 1)) %>%
  mutate(total_response_count_lag1 = dplyr::lag(total_response_count, n = 1)) %>%
  mutate(ratio_lag1 = dplyr::lag(ratio, n = 1)) %>%
  mutate(avg_mip_country_immigration_lag1 = dplyr::lag(avg_mip_country_immigration, n = 1)) %>%
  mutate(weighted_avg_mip_country_immigration_lag1 = dplyr::lag(weighted_avg_mip_country_immigration, n = 1)) %>%
  mutate(monthly_resp_f_count_lag2 = dplyr::lag(monthly_resp_f_count, n = 2)) %>%
  mutate(monthly_tp_h_count_lag2 = dplyr::lag(monthly_tp_h_count, n = 2)) %>%
  mutate(asylum_total_lag2 = dplyr::lag(asylum_total, n = 2)) %>%
  mutate(asylum_total_scaled_lag2 = dplyr::lag(asylum_total_scaled, n = 2)) %>%
  mutate(total_response_count_lag2 = dplyr::lag(total_response_count, n = 2)) %>%
  mutate(ratio_lag2 = dplyr::lag(ratio, n = 2)) %>%
  mutate(avg_mip_country_immigration_lag2 = dplyr::lag(avg_mip_country_immigration, n = 2)) %>%
  mutate(weighted_avg_mip_country_immigration_lag2 = dplyr::lag(weighted_avg_mip_country_immigration, n = 2)) %>%
  mutate(monthly_resp_f_count_lag3 = dplyr::lag(monthly_resp_f_count, n = 3)) %>%
  mutate(monthly_tp_h_count_lag3 = dplyr::lag(monthly_tp_h_count, n = 3)) %>%
  mutate(asylum_total_lag3 = dplyr::lag(asylum_total, n = 3)) %>%
  mutate(asylum_total_scaled_lag3 = dplyr::lag(asylum_total_scaled, n = 3)) %>%
  mutate(total_response_count_lag3 = dplyr::lag(total_response_count, n = 3)) %>%
  mutate(ratio_lag3 = dplyr::lag(ratio, n = 3)) %>%
  mutate(avg_mip_country_immigration_lag3 = dplyr::lag(avg_mip_country_immigration, n = 3)) %>%
  mutate(weighted_avg_mip_country_immigration_lag3 = dplyr::lag(weighted_avg_mip_country_immigration, n = 3)) %>%
  mutate(monthly_resp_f_count_lag6 = dplyr::lag(monthly_resp_f_count, n = 6)) %>%
  mutate(monthly_tp_h_count_lag6 = dplyr::lag(monthly_tp_h_count, n = 6)) %>%
  mutate(asylum_total_lag6 = dplyr::lag(asylum_total, n = 6)) %>%
  mutate(asylum_total_scaled_lag6 = dplyr::lag(asylum_total_scaled, n = 6)) %>%
  mutate(total_response_count_lag6 = dplyr::lag(total_response_count, n = 6)) %>%
  mutate(ratio_lag6 = dplyr::lag(ratio, n = 6)) %>%
  mutate(avg_mip_country_immigration_lag6 = dplyr::lag(avg_mip_country_immigration, n = 6)) %>%
  mutate(weighted_avg_mip_country_immigration_lag6 = dplyr::lag(weighted_avg_mip_country_immigration, n = 6)) %>%
  mutate(monthly_resp_f_count_lag12 = dplyr::lag(monthly_resp_f_count, n = 12)) %>%
  mutate(monthly_tp_h_count_lag12 = dplyr::lag(monthly_tp_h_count, n = 12)) %>%
  mutate(asylum_total_lag12 = dplyr::lag(asylum_total, n = 12)) %>%
  mutate(asylum_total_scaled_lag12 = dplyr::lag(asylum_total_scaled, n = 12)) %>%
  mutate(total_response_count_lag12 = dplyr::lag(total_response_count, n = 12)) %>%
  mutate(ratio_lag12 = dplyr::lag(ratio, n = 12)) %>%
  mutate(avg_mip_country_immigration_lag12 = dplyr::lag(avg_mip_country_immigration, n = 12)) %>%
  mutate(weighted_avg_mip_country_immigration_lag12 = dplyr::lag(weighted_avg_mip_country_immigration, n = 12))

# add post-Lisbon variable
main_data$post_lisbon <- ifelse(main_data$month >= "2009-12-01", 1, 0)

```

\newpage

# Descriptives

```{r, descriptives, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

df_descriptive <- pivot_longer(main_data %>%
                    mutate(asylum_total_new = asylum_total / 10000) %>%
                    mutate(avg_mip_country_immigration_new = avg_mip_country_immigration * 100) %>%
                    mutate(
                      total_response_count_avg = rollapply(
                        total_response_count,
                        width = 12, FUN = mean, align = "right", fill = NA)
                      ),
                    cols = c("total_response_count_avg",
                             "asylum_total_new",
                             "avg_mip_country_immigration_new"),
                    names_to = "variable") %>%
  dplyr::select("month", "variable", "value") %>%
  mutate(month = as.Date(month)) %>%
  mutate(variable = recode(variable, 
                           "total_response_count_avg" = "Time pressure invocations (12-month average)",
                           "asylum_total_new" = "Asylum numbers",
                           "avg_mip_country_immigration_new" = "Public pressure"))

shaded_areas <- data.frame(
  start = as.Date(c('1993-11-01', '1999-05-01', '2009-12-01')),
  end = as.Date(c('1999-04-30', '2009-11-30', '2023-12-31')),
  color = c('lightblue', 'lightgreen', 'lightcoral'),
  label = c('Maastricht', 'Amsterdam', 'Lisbon')
)

# events <- data.frame(
#   year = as.Date(c('1995-01-01', '1999-01-01', '2004-01-01', 
#                    '2008-01-01', '2010-01-01', '2014-01-01', '2019-01-01')),
#   label = c('Gradin', 'Vitorino', 'Fratini', 'Barrot', 
#             'Malmström', 'Avramopoulos', 'Johansson')
# )

ggplot(df_descriptive, aes(x = month, y = value, colour = variable)) +
  geom_point(alpha = 1) +
  geom_smooth() + #smoothed line
  #geom_line() + #line plot
  geom_rect(data = shaded_areas, 
            aes(xmin = start, xmax = end, ymin = -Inf, ymax = Inf, fill = color), 
            alpha = 0.25, inherit.aes = FALSE) +
  annotate("text", 
           x = as.Date(c('1993-11-01', '1999-05-01', '2009-12-01')),  # Start + 1 month
           y = 40,
           label = shaded_areas$label, 
           size = 4, color = 'black', hjust = -0.1, vjust = 1) +
  geom_vline(xintercept = as.Date(c('1993-11-01', '1999-05-01', '2009-12-01')), linetype = "dashed") +
  scale_color_manual(values = c("#000099", "#FFCC00", "#000000")) +
  scale_fill_identity() +
  scale_x_date(limits = c(as.Date("1990-01-01"), as.Date("2023-12-31")), expand = c(0,5)) +
  scale_y_continuous(limits = c(0, 45), expand = c(0,0)) +  # Ensure y-axis starts at 0
  labs(x = "", y = "", colour = "",
       caption = "Note: The points represent the monthly values for asylum numbers, public pressure and time pressure \ninvocations. The lines represent the smoothed values for the respective variables. In the case of time \npressure invocations, the 12-month moving average values are reported.") +
  theme_bw() +
  theme(legend.position = "top",
        plot.caption = element_text(face = "italic", hjust=0))

```

\newpage

# VAR

## Models

```{r, var - time series data, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

ts_response <- ts(
  main_data$total_response_count,
  start = c(1990, 1), end = c(2023, 12), frequency = 12
  )

#PP.test(ts_response)

ts_asylum <- ts(
  main_data$asylum_total_scaled,
  start = c(1990, 1), end = c(2023, 12), frequency = 12
  )

#PP.test(ts_asylum)

ts_mip <- ts(
  main_data$avg_mip_country_immigration,
  start = c(1990, 1), end = c(2023, 12), frequency = 12
  )

#PP.test(ts_mip)

post_lisbon_full <- main_data %>%
  dplyr::select(post_lisbon) %>%
  unlist() %>%
  as.numeric()

post_lisbon_short <- main_data %>%
  filter(month >= "2002-06-01" & month < "2024-01-01") %>% 
  dplyr::select(post_lisbon) %>%
  unlist() %>%
  as.numeric()

```

```{r, var - models, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# full time period

var_data_1 <- window(
  ts.union(ts_response, ts_asylum), 
  start = c(1990, 1), end = c(2023, 12)
  )
var_1 <- VAR(scale(var_data_1), type = "const", exogen = post_lisbon_full,
             lag.max = 12, ic = "AIC")

var_data_2 <- window(
  ts.union(ts_response, ts_mip), 
  start = c(2002, 6), end = c(2023, 12)
  )
var_2 <- VAR(scale(var_data_2), type = "const", exogen = post_lisbon_short,
             lag.max = 12, ic = "AIC")

var_data_3 <- window(
  ts.union(ts_asylum, ts_mip), 
  start = c(2002, 6), end = c(2023, 12)
  )
var_3 <- VAR(scale(var_data_3), type = "const", exogen = post_lisbon_short,
             lag.max = 12, ic = "AIC")

var_data_4 <- window(
  ts.union(ts_response, ts_asylum, ts_mip), 
  start = c(2002, 6), end = c(2023, 12)
  )
var_4 <- VAR(scale(var_data_4), type = "const", exogen = post_lisbon_short,
             lag.max = 12, ic = "AIC")

# pre-lisbon period

var_data_5 <- window(
  ts.union(ts_response, ts_asylum), 
  start = c(1990, 01), end = c(2009, 11)
  )
var_5 <- VAR(scale(var_data_5), type = "const",
             lag.max = 12, ic = "AIC")

var_data_6 <- window(
  ts.union(ts_response, ts_mip), 
  start = c(2002, 06), end = c(2009, 11)
  )
var_6 <- VAR(scale(var_data_6), type = "const",
             lag.max = 12, ic = "AIC")

var_data_7 <- window(
  ts.union(ts_asylum, ts_mip), 
  start = c(2002, 06), end = c(2009, 11)
  )
var_7 <- VAR(scale(var_data_7), type = "const",
             lag.max = 12, ic = "AIC")

var_data_8 <- window(
  ts.union(ts_response, ts_asylum, ts_mip), 
  start = c(2002, 6), end = c(2009, 11)
  )
var_8 <- VAR(scale(var_data_8), type = "const",
             lag.max = 12, ic = "AIC")

# post-lisbon period

var_data_9 <- window(
  ts.union(ts_response, ts_asylum), 
  start = c(2009, 12), end = c(2023, 12)
  )
var_9 <- VAR(scale(var_data_9), type = "const",
             lag.max = 12, ic = "AIC")

var_data_10 <- window(
  ts.union(ts_response, ts_mip), 
  start = c(2009, 12), end = c(2023, 12)
  )
var_10 <- VAR(scale(var_data_10), type = "const",
              lag.max = 12, ic = "AIC")

var_data_11 <- window(
  ts.union(ts_asylum, ts_mip), 
  start = c(2009, 12), end = c(2023, 12)
  )
var_11 <- VAR(scale(var_data_11), type = "const",
              lag.max = 12, ic = "AIC")

var_data_12 <- window(
  ts.union(ts_response, ts_asylum, ts_mip), 
  start = c(2009, 12), end = c(2023, 12)
  )
var_12 <- VAR(scale(var_data_12), type = "const",
              lag.max = 12, ic = "AIC")

```

\newpage

## 12-month cumulative IRF

```{r, var - irf, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# full time period

## response: time pressure invocations; impulse: asylum numbers
irf_1 <- irf(var_1, impulse="ts_asylum", response="ts_response", ortho = FALSE,
             cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_2 <- irf(var_4, impulse="ts_asylum", response="ts_response", ortho = FALSE,
             cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: time pressure invocations; impulse: mip
irf_3 <- irf(var_2, impulse="ts_mip", response="ts_response", ortho = FALSE,
             cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_4 <- irf(var_4, impulse="ts_mip", response="ts_response", ortho = FALSE,
             cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: asylum numbers; impulse: time pressure invocations
irf_5 <- irf(var_1, impulse="ts_response", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # bivariate model
irf_6 <- irf(var_4, impulse="ts_response", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # multivariate model

## response: asylum numbers; impulse: mip
irf_7 <- irf(var_3, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # bivariate model
irf_8 <- irf(var_4, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: time pressure invocations
irf_9 <- irf(var_2, impulse="ts_response", response="ts_mip", ortho = FALSE,
             cumulative = TRUE, seed = 1) # bivariate model
irf_10 <- irf(var_4, impulse="ts_response", response="ts_mip", ortho = FALSE,
              cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: asylum numbers
irf_11 <- irf(var_3, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
              cumulative = TRUE, seed = 1) # bivariate model
irf_12 <- irf(var_4, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
              cumulative = TRUE, seed = 1) # multivariate model

# pre-lisbon period

## response: time pressure invocations; impulse: asylum numbers
irf_13 <- irf(var_5, impulse="ts_asylum", response="ts_response", ortho = FALSE,
              cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_14 <- irf(var_8, impulse="ts_asylum", response="ts_response", ortho = FALSE,
              cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: time pressure invocations; impulse: mip
irf_15 <- irf(var_6, impulse="ts_mip", response="ts_response", ortho = FALSE,
             cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_16 <- irf(var_8, impulse="ts_mip", response="ts_response", ortho = FALSE,
             cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: asylum numbers; impulse: time pressure invocations
irf_17 <- irf(var_5, impulse="ts_response", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # bivariate model
irf_18 <- irf(var_8, impulse="ts_response", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # multivariate model

## response: asylum numbers; impulse: mip
irf_19 <- irf(var_8, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # bivariate model
irf_20 <- irf(var_8, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: time pressure invocations
irf_21 <- irf(var_6, impulse="ts_response", response="ts_mip", ortho = FALSE,
             cumulative = TRUE, seed = 1) # bivariate model
irf_22 <- irf(var_8, impulse="ts_response", response="ts_mip", ortho = FALSE,
              cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: asylum numbers
irf_23 <- irf(var_7, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
              cumulative = TRUE, seed = 1) # bivariate model
irf_24 <- irf(var_8, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
              cumulative = TRUE, seed = 1) # multivariate model

# post-lisbon period

## response: time pressure invocations; impulse: asylum numbers
irf_25 <- irf(var_9, impulse="ts_asylum", response="ts_response", ortho = FALSE,
              cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_26 <- irf(var_12, impulse="ts_asylum", response="ts_response", ortho = FALSE,
              cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: time pressure invocations; impulse: mip
irf_27 <- irf(var_10, impulse="ts_mip", response="ts_response", ortho = FALSE,
             cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_28 <- irf(var_12, impulse="ts_mip", response="ts_response", ortho = FALSE,
             cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: asylum numbers; impulse: time pressure invocations
irf_29 <- irf(var_9, impulse="ts_response", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # bivariate model
irf_30 <- irf(var_12, impulse="ts_response", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # multivariate model

## response: asylum numbers; impulse: mip
irf_31 <- irf(var_11, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # bivariate model
irf_32 <- irf(var_12, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
             cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: time pressure invocations
irf_33 <- irf(var_10, impulse="ts_response", response="ts_mip", ortho = FALSE,
             cumulative = TRUE, seed = 1) # bivariate model
irf_34 <- irf(var_12, impulse="ts_response", response="ts_mip", ortho = FALSE,
              cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: asylum numbers
irf_35 <- irf(var_11, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
              cumulative = TRUE, seed = 1) # bivariate model
irf_36 <- irf(var_12, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
              cumulative = TRUE, seed = 1) # multivariate model

```

\newpage

```{r, var - irf graphics, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# data frame
irf_df <- data.frame(
  var_model = c("bivariate", "multivariate", "bivariate", "multivariate", 
                "bivariate", "multivariate", "bivariate", "multivariate",
                "bivariate", "multivariate", "bivariate", "multivariate",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)"),
  response = c(irf_1$response, irf_2$response, irf_3$response, irf_4$response,
               irf_5$response, irf_6$response, irf_7$response, irf_8$response,
               irf_9$response, irf_10$response, irf_11$response, irf_12$response,
               irf_13$response, irf_14$response, irf_15$response, irf_16$response,
               irf_17$response, irf_18$response, irf_19$response, irf_20$response,
               irf_21$response, irf_22$response, irf_23$response, irf_24$response,
               irf_25$response, irf_26$response, irf_27$response, irf_28$response,
               irf_29$response, irf_30$response, irf_31$response, irf_32$response,
               irf_33$response, irf_34$response, irf_35$response, irf_36$response),
  impulse = c(irf_1$impulse, irf_2$impulse, irf_3$impulse, irf_4$impulse,
              irf_5$impulse, irf_6$impulse, irf_7$impulse, irf_8$impulse,
              irf_9$impulse, irf_10$impulse, irf_11$impulse, irf_12$impulse,
              irf_13$impulse, irf_14$impulse, irf_15$impulse, irf_16$impulse,
              irf_17$impulse, irf_18$impulse, irf_19$impulse, irf_20$impulse,
              irf_21$impulse, irf_22$impulse, irf_23$impulse, irf_24$impulse,
              irf_25$impulse, irf_26$impulse, irf_27$impulse, irf_28$impulse,
              irf_29$impulse, irf_30$impulse, irf_31$impulse, irf_32$impulse,
              irf_33$impulse, irf_34$impulse, irf_35$impulse, irf_36$impulse),
  irf = c(tail(unlist(irf_1$irf), n = 1), tail(unlist(irf_2$irf), n = 1),
          tail(unlist(irf_3$irf), n = 1), tail(unlist(irf_4$irf), n = 1),
          tail(unlist(irf_5$irf), n = 1), tail(unlist(irf_6$irf), n = 1),
          tail(unlist(irf_7$irf), n = 1), tail(unlist(irf_8$irf), n = 1),
          tail(unlist(irf_9$irf), n = 1), tail(unlist(irf_10$irf), n = 1),
          tail(unlist(irf_11$irf), n = 1), tail(unlist(irf_12$irf), n = 1),
          tail(unlist(irf_13$irf), n = 1), tail(unlist(irf_14$irf), n = 1),
          tail(unlist(irf_15$irf), n = 1), tail(unlist(irf_16$irf), n = 1),
          tail(unlist(irf_17$irf), n = 1), tail(unlist(irf_18$irf), n = 1),
          tail(unlist(irf_19$irf), n = 1), tail(unlist(irf_20$irf), n = 1),
          tail(unlist(irf_21$irf), n = 1), tail(unlist(irf_22$irf), n = 1),
          tail(unlist(irf_23$irf), n = 1), tail(unlist(irf_24$irf), n = 1),
          tail(unlist(irf_25$irf), n = 1), tail(unlist(irf_26$irf), n = 1),
          tail(unlist(irf_27$irf), n = 1), tail(unlist(irf_28$irf), n = 1),
          tail(unlist(irf_29$irf), n = 1), tail(unlist(irf_30$irf), n = 1),
          tail(unlist(irf_31$irf), n = 1), tail(unlist(irf_32$irf), n = 1),
          tail(unlist(irf_33$irf), n = 1), tail(unlist(irf_34$irf), n = 1),
          tail(unlist(irf_35$irf), n = 1), tail(unlist(irf_36$irf), n = 1)),
  lower = c(tail(unlist(irf_1$Lower), n = 1), tail(unlist(irf_2$Lower), n = 1),
            tail(unlist(irf_3$Lower), n = 1), tail(unlist(irf_4$Lower), n = 1),
            tail(unlist(irf_5$Lower), n = 1), tail(unlist(irf_6$Lower), n = 1),
            tail(unlist(irf_7$Lower), n = 1), tail(unlist(irf_8$Lower), n = 1),
            tail(unlist(irf_9$Lower), n = 1), tail(unlist(irf_10$Lower), n = 1),
            tail(unlist(irf_11$Lower), n = 1), tail(unlist(irf_12$Lower), n = 1),
            tail(unlist(irf_13$Lower), n = 1), tail(unlist(irf_14$Lower), n = 1),
            tail(unlist(irf_15$Lower), n = 1), tail(unlist(irf_16$Lower), n = 1),
            tail(unlist(irf_17$Lower), n = 1), tail(unlist(irf_18$Lower), n = 1),
            tail(unlist(irf_19$Lower), n = 1), tail(unlist(irf_20$Lower), n = 1),
            tail(unlist(irf_21$Lower), n = 1), tail(unlist(irf_22$Lower), n = 1),
            tail(unlist(irf_23$Lower), n = 1), tail(unlist(irf_24$Lower), n = 1),
            tail(unlist(irf_25$Lower), n = 1), tail(unlist(irf_26$Lower), n = 1),
            tail(unlist(irf_27$Lower), n = 1), tail(unlist(irf_28$Lower), n = 1),
            tail(unlist(irf_29$Lower), n = 1), tail(unlist(irf_30$Lower), n = 1),
            tail(unlist(irf_31$Lower), n = 1), tail(unlist(irf_32$Lower), n = 1),
            tail(unlist(irf_33$Lower), n = 1), tail(unlist(irf_34$Lower), n = 1),
            tail(unlist(irf_35$Lower), n = 1), tail(unlist(irf_36$Lower), n = 1)),
  upper = c(tail(unlist(irf_1$Upper), n = 1), tail(unlist(irf_2$Upper), n = 1),
            tail(unlist(irf_3$Upper), n = 1), tail(unlist(irf_4$Upper), n = 1),
            tail(unlist(irf_5$Upper), n = 1), tail(unlist(irf_6$Upper), n = 1),
            tail(unlist(irf_7$Upper), n = 1), tail(unlist(irf_8$Upper), n = 1),
            tail(unlist(irf_9$Upper), n = 1), tail(unlist(irf_10$Upper), n = 1),
            tail(unlist(irf_11$Upper), n = 1), tail(unlist(irf_12$Upper), n = 1),
            tail(unlist(irf_13$Upper), n = 1), tail(unlist(irf_14$Upper), n = 1),
            tail(unlist(irf_15$Upper), n = 1), tail(unlist(irf_16$Upper), n = 1),
            tail(unlist(irf_17$Upper), n = 1), tail(unlist(irf_18$Upper), n = 1),
            tail(unlist(irf_19$Upper), n = 1), tail(unlist(irf_20$Upper), n = 1),
            tail(unlist(irf_21$Upper), n = 1), tail(unlist(irf_22$Upper), n = 1),
            tail(unlist(irf_23$Upper), n = 1), tail(unlist(irf_24$Upper), n = 1),
            tail(unlist(irf_25$Upper), n = 1), tail(unlist(irf_26$Upper), n = 1),
            tail(unlist(irf_27$Upper), n = 1), tail(unlist(irf_28$Upper), n = 1),
            tail(unlist(irf_29$Upper), n = 1), tail(unlist(irf_30$Upper), n = 1),
            tail(unlist(irf_31$Upper), n = 1), tail(unlist(irf_32$Upper), n = 1),
            tail(unlist(irf_33$Upper), n = 1), tail(unlist(irf_34$Upper), n = 1),
            tail(unlist(irf_35$Upper), n = 1), tail(unlist(irf_36$Upper), n = 1))
  ) %>%
  mutate(response = recode(response, 
                           "ts_response" = "Time pressure invocations",
                           "ts_asylum" = "Asylum numbers",
                           "ts_mip" = "Public pressure")) %>%
  mutate(impulse = recode(impulse, 
                          "ts_response" = "Time pressure invocations",
                          "ts_asylum" = "Asylum numbers",
                          "ts_mip" = "Public pressure")) %>%
  mutate(var_model =  fct_relevel(var_model,
                                  "bivariate", "multivariate",
                                  "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                                  "bivariate (post-lisbon)", "multivariate (post-lisbon)"))

# plots
ggplot(irf_df %>% filter(response == "Time pressure invocations"), 
       aes(x = fct_rev(impulse), y = irf, ymin = lower, ymax = upper, colour = var_model)) +
  geom_linerange(size = 2.5, position = position_dodge(width = -0.55), alpha = 0.4) +  
  geom_point(position = position_dodge(width = -0.55), size = 3.5) +
  geom_hline(aes(yintercept = 0), color = "black", linetype = 2, alpha = 0.75) +
  scale_color_manual(values=c("#000000", "#999999", "#000099", "#0099FF", "#FFCC00", "#FFFF33")) +
  coord_flip() +
  #scale_y_continuous(limits = c(-2, 6), breaks = c(-2, 0, 2, 4, 6)) +
  labs(subtitle = "",
       title = "",
       y = "Predicted time pressure invocations (12 months cumulative)", x = "",
       colour = "VAR Model") +
  theme_bw()

ggplot(irf_df %>% filter(response == "Asylum numbers"), 
       aes(x = fct_rev(impulse), y = irf, ymin = lower, ymax = upper, colour = var_model)) +
  geom_linerange(size = 2.5, position = position_dodge(width = -0.55), alpha = 0.4) +  
  geom_point(position = position_dodge(width = -0.55), size = 3.5) +
  geom_hline(aes(yintercept = 0), color = "black", linetype = 2, alpha = 0.75) +
  scale_x_discrete(
    labels = c("Time pressure\ninvocations", "Public pressure")
  ) +
  scale_color_manual(values=c("#000000", "#999999", "#000099", "#0099FF", "#FFCC00", "#FFFF33")) +
  coord_flip() +
  #scale_y_continuous(limits = c(-2, 6), breaks = c(-2, 0, 2, 4, 6)) +
  labs(subtitle = "",
       title = "",
       y = "Predicted asylum numbers (12 months cumulative)", x = "",
       colour = "VAR Model") +
  theme_bw()

ggplot(irf_df %>% filter(response == "Public pressure"), 
       aes(x = fct_rev(impulse), y = irf, ymin = lower, ymax = upper, colour = var_model)) +
  geom_linerange(size = 2.5, position = position_dodge(width = -0.55), alpha = 0.4) +  
  geom_point(position = position_dodge(width = -0.55), size = 3.5) +
  geom_hline(aes(yintercept = 0), color = "black", linetype = 2, alpha = 0.75) +
  scale_x_discrete(
    labels = c("Time pressure\ninvocations", "Asylum numbers")
  ) +
  scale_color_manual(values=c("#000000", "#999999", "#000099", "#0099FF", "#FFCC00", "#FFFF33")) +
  coord_flip() +
  #scale_y_continuous(limits = c(-2, 6), breaks = c(-2, 0, 2, 4, 6)) +
  labs(subtitle = "",
       title = "",
       y = "Predicted public pressure (12 months cumulative)", x = "",
       colour = "VAR Model") +
  theme_bw()

```

\newpage

# Mediation Analysis

```{r}

# Install packages if necessary
install.packages("mediation")
library(mediation)

```

```{r}
# Fit the mediator model (Mediator = weighted_av_mip_country_immigration)
mediator_model <- lm(weighted_avg_mip_country_immigration ~ asylum_total_scaled, data = main_data)

# Fit the outcome model (DV = total_response_count)
outcome_model <- lm(total_response_count ~ asylum_total_scaled + weighted_avg_mip_country_immigration, data = main_data)


```

\# Perform mediation analysis

mediation_result \<- mediate(mediator_model, outcome_model,

treat = "asylum_total_scaled", mediator = "weighted_avg_mip_country_immigration",

boot = TRUE, sims = 500) \# Bootstrapping for confidence intervals

\# View the results

summary(mediation_result)

```{r}
# Perform mediation analysis
mediation_result <- mediate(mediator_model, outcome_model, 
                            treat = "asylum_total_scaled", mediator = "weighted_avg_mip_country_immigration", 
                            boot = TRUE, sims = 500)  # Bootstrapping for confidence intervals

# View the results
summary(mediation_result)

```

# Robustness checks - MY confounding

```{r}
# Perform sensitivity analysis
sensitivity_result <- medsens(mediation_result)

# View the sensitivity results
summary(sensitivity_result)

```

## Mediation for post-Lisbon

```{r}
# Extract the year from the 'month' column
main_data$year <- format(main_data$month, "%Y")
main_data$year <- as.numeric(main_data$year)  # Convert it to numeric


```

```{r}
# Filter the data to include only years greater than 2009
filtered_data <- main_data[main_data$year > 2009, ]

# Fit the mediator model (mediator as a function of treatment) using the filtered data
mediator_model <- lm(weighted_avg_mip_country_immigration ~ asylum_total_scaled, data = filtered_data)

# Fit the outcome model (outcome as a function of treatment and mediator) using the filtered data
outcome_model <- lm(total_response_count ~ asylum_total_scaled + weighted_avg_mip_country_immigration, data = filtered_data)

# Perform mediation analysis with the filtered data
mediation_result <- mediate(mediator_model, outcome_model, 
                            treat = "asylum_total_scaled", mediator = "weighted_avg_mip_country_immigration", 
                            boot = TRUE, sims = 500)  # Bootstrapping for confidence intervals

# View the results
summary(mediation_result)

```

```{r}
# Perform sensitivity analysis
sensitivity_result <- medsens(mediation_result)

# View the sensitivity results
summary(sensitivity_result)
```

\newpage

# Appendix

## A.1 Public opinion data: Descriptives for "Most important problem" in Eurobarometer data

### Eurobarometer surveys

To construct the independent variable for public pressure, we rely on citizens' responses to the questions on the most important problem (MIP) that their country is currently facing in Eurobarometer surveys (2002-2023). The following table provides an overview of the Eurobarometer surveys, which contain the MIP item and were used to construct the independent variable.

```{r, eurobarometer data, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# load data
load("../Data/eurobarometer_mip_immigration_data.RDA")

mip_var_report %>%
  dplyr::select(studyno, eurobarometer, fieldwork_month, fieldwork_year) %>%
  knitr::kable(col.names = c("Study number", "Eurobarometer", 
                             "Fieldwork month(s)", "Fieldwork year"),
               row.names = TRUE,
               caption = "Overview of Eurobarometer surveys used to construct the
               independent variable for public pressure.",
               align = c("c", "c", "c", "c"))

```

### Descriptives

As discussed in the article and above, the independent variable for public pressure is based on citizens' responses to the question on the 'most important problem' item in Eurobarometer surveys. In the article, we rely on the item related the most important problem (MIP) facing the respective country of a respondent. More precisely, we calculate the average importance of immigration in response to the MIP question per survey. The standardised MIP item regularly exists in Eurobarometer surveys since 2002 and is available at a relatively high granularity (typically two to three surveys contain the MIP item per year). Since 2010 Eurobarometer surveys also regularly contain a standardised item related to the most important problem (MIP) faced by the European Union. However, in the article, we opted to use the MIP for country to construct the independent variable for public pressure because it allows to investigate a longer time period at a comparatively high granularity.

The following figures provide descriptives on the development of public pressure related to migration and asylum based on the Eurobarometer data. In addtion, the figures compare the responses to the MIP items for country and for the EU as well as weighted and unweighted values. Overall, the development is largely parallel for all measures. However, the values related to MIP for the EU are in general higher than in the case of MIP for country, particularly during the migration crisis in 2015 and 2016.

```{r, mip descriptives, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# salience (without weighting)
overview_salience_unweighted <- left_join(
  
  left_join(
    mip_data %>%
      filter(!is.na(weight_EU)) %>%
      filter(weight_EU > 0) %>%
      filter(!is.na(mip_country_immigration)) %>%
      group_by(fieldwork_year) %>%
      summarise(sum_country_unweighted = sum(mip_country_immigration)),
    mip_data %>%
      filter(!is.na(weight_EU)) %>%
      filter(weight_EU > 0) %>%
      filter(!is.na(mip_country_immigration)) %>%
      group_by(fieldwork_year) %>%
      count()
  ) %>%
    mutate(salience_country_unweighted = sum_country_unweighted / n * 100) %>%
    dplyr::select(fieldwork_year, salience_country_unweighted),
  
  left_join(
    mip_data %>%
      filter(!is.na(weight_EU)) %>%
      filter(weight_EU > 0) %>%
      filter(!is.na(mip_eu_immigration)) %>%
      group_by(fieldwork_year) %>%
      summarise(sum_eu_unweighted = sum(mip_eu_immigration)),
    mip_data %>%
      filter(!is.na(weight_EU)) %>%
      filter(weight_EU > 0) %>%
      filter(!is.na(mip_eu_immigration)) %>%
      group_by(fieldwork_year) %>%
      count()
  ) %>%
    mutate(salience_eu_unweighted = sum_eu_unweighted / n * 100) %>%
    dplyr::select(fieldwork_year, salience_eu_unweighted)
  
)

plot_unweighted <- overview_salience_unweighted %>%
  pivot_longer(cols = starts_with("salience_"),
               names_to = "issue",
               names_prefix = "salience_",
               values_to = "salience") %>%
  ggplot(aes(x = fieldwork_year, y = salience, color = issue)) +
  geom_point() +
  geom_line() + 
  scale_x_continuous(limits = c(2002, 2023)) +
  scale_y_continuous(limits = c(0, 52)) +
  labs(x = "Fieldwork year", y = "Public pressure \n(unweighted)") +
  scale_color_manual(name="Level", 
                     labels=c("Country", "EU"), 
                     values=c("red", "blue")) +
  theme_bw()


## salience (with weighting)
overview_salience_weighted <- left_join(
  
  left_join(
    mip_data %>%
      filter(!is.na(weight_EU)) %>%
      filter(weight_EU > 0) %>%
      filter(!is.na(mip_country_immigration)) %>%
      mutate(mip_country_weighted = mip_country_immigration * weight_EU) %>%
      group_by(fieldwork_year) %>%
      summarise(sum_country_weighted = sum(mip_country_weighted)),
    mip_data %>%
      filter(!is.na(weight_EU)) %>%
      filter(weight_EU > 0) %>%
      filter(!is.na(mip_country_immigration)) %>%
      group_by(fieldwork_year) %>%
      count()
  ) %>%
    mutate(salience_country_weighted = sum_country_weighted / n * 100) %>%
    dplyr::select(fieldwork_year, salience_country_weighted),
  
  left_join(
    mip_data %>%
      filter(!is.na(weight_EU)) %>%
      filter(weight_EU > 0) %>%
      filter(!is.na(mip_eu_immigration)) %>%
      mutate(mip_eu_weighted = mip_eu_immigration * weight_EU) %>%
      group_by(fieldwork_year) %>%
      summarise(sum_eu_weighted = sum(mip_eu_weighted)),
    mip_data %>%
      filter(!is.na(weight_EU)) %>%
      filter(weight_EU > 0) %>%
      filter(!is.na(mip_eu_immigration)) %>%
      group_by(fieldwork_year) %>%
      count()
  ) %>%
    mutate(salience_eu_weighted = sum_eu_weighted / n * 100) %>%
    dplyr::select(fieldwork_year, salience_eu_weighted)
  
)

plot_weighted <- overview_salience_weighted %>%
  pivot_longer(cols = starts_with("salience_"),
               names_to = "issue",
               names_prefix = "salience_",
               values_to = "salience") %>%
  ggplot(aes(x = fieldwork_year, y = salience, color = issue)) +
  geom_point() +
  geom_line() + 
  scale_x_continuous(limits = c(2002, 2023)) +
  scale_y_continuous(limits = c(0, 52)) +
  labs(x = "Fieldwork year", y = "Public pressure \n(weighted)") +
  scale_color_manual(name="Level", 
                     labels=c("Country", "EU"), 
                     values=c("red", "blue")) +
  theme_bw()

## plot
grid.arrange(plot_unweighted, plot_weighted)

```

\newpage

## A.2 VECM: Vector Error Correction Models

Vector Error Corrections models (VECMs) incorporate error correction terms and are useful for handling non-stationary time series. Both features allow to control or correct for common problems faced by Vector Autoregression (VAR) models, such as stationary/non-stationary time series and deviations from long-term equilibrium. Therefore, we replicate the analysis reported in the article by using VECMs. The main results from this robustness check are largely consistent with the findings reported in the article.

### Models

```{r, vec - models, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# full time period
vec_a1 <- ca.jo(scale(var_data_1), ecdet = "none", type = "trace",
               K = 12, spec = "transitory", dumvar = post_lisbon_full)
var_a1 <- vec2var(vec_a1, r = 1)

vec_a2 <- ca.jo(scale(var_data_2), ecdet = "none", type = "trace",
               K = 6, spec = "transitory", dumvar = post_lisbon_short)
var_a2 <- vec2var(vec_a2, r = 1)

vec_a3 <- ca.jo(scale(var_data_3), ecdet = "none", type = "trace",
               K = 6, spec = "transitory", dumvar = post_lisbon_short)
var_a3 <- vec2var(vec_a3, r = 1)

vec_a4 <- ca.jo(scale(var_data_4), ecdet = "none", type = "trace",
               K = 6, spec = "transitory", dumvar = post_lisbon_short)
var_a4 <- vec2var(vec_a4, r = 2)

# pre-lisbon period
vec_a5 <- ca.jo(scale(var_data_5), ecdet = "none", type = "trace",
               K = 12, spec = "transitory")
var_a5 <- vec2var(vec_a5, r = 1)

vec_a6 <- ca.jo(scale(var_data_6), ecdet = "none", type = "trace",
               K = 6, spec = "transitory")
var_a6 <- vec2var(vec_a6, r = 1)

vec_a7 <- ca.jo(scale(var_data_7), ecdet = "none", type = "trace",
               K = 6, spec = "transitory")
var_a7 <- vec2var(vec_a7, r = 1)

vec_a8 <- ca.jo(scale(var_data_8), ecdet = "none", type = "trace",
               K = 6, spec = "transitory")
var_a8 <- vec2var(vec_a8, r = 2)

# post-lisbon period
vec_a9 <- ca.jo(scale(var_data_9), ecdet = "none", type = "trace",
               K = 12, spec = "transitory")
var_a9 <- vec2var(vec_a9, r = 1)

vec_a10 <- ca.jo(scale(var_data_10), ecdet = "none", type = "trace",
               K = 6, spec = "transitory")
var_a10 <- vec2var(vec_a10, r = 1)

vec_a11 <- ca.jo(scale(var_data_11), ecdet = "none", type = "trace",
               K = 6, spec = "transitory")
var_a11 <- vec2var(vec_a11, r = 1)

vec_a12 <- ca.jo(scale(var_data_12), ecdet = "none", type = "trace",
               K = 6, spec = "transitory")
var_a12 <- vec2var(vec_a12, r = 2)

```

\newpage

### 12-month cumulative IRF

```{r, vec - irf, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# full time period

## response: time pressure invocations; impulse: asylum numbers
irf_a1 <- irf(var_a1, impulse="ts_asylum", response="ts_response", ortho = FALSE,
              cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_a2 <- irf(var_a4, impulse="ts_asylum", response="ts_response", ortho = FALSE,
              cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: time pressure invocations; impulse: mip
irf_a3 <- irf(var_a2, impulse="ts_mip", response="ts_response", ortho = FALSE,
              cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_a4 <- irf(var_a4, impulse="ts_mip", response="ts_response", ortho = FALSE,
              cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: asylum numbers; impulse: time pressure invocations
irf_a5 <- irf(var_a1, impulse="ts_response", response="ts_asylum", ortho = FALSE,
              cumulative = TRUE, seed = 1) # bivariate model
irf_a6 <- irf(var_a4, impulse="ts_response", response="ts_asylum", ortho = FALSE,
              cumulative = TRUE, seed = 1) # multivariate model

## response: asylum numbers; impulse: mip
irf_a7 <- irf(var_a3, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
              cumulative = TRUE, seed = 1) # bivariate model
irf_a8 <- irf(var_a4, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
              cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: time pressure invocations
irf_a9 <- irf(var_a2, impulse="ts_response", response="ts_mip", ortho = FALSE,
              cumulative = TRUE, seed = 1) # bivariate model
irf_a10 <- irf(var_a4, impulse="ts_response", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: asylum numbers
irf_a11 <- irf(var_a3, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # bivariate model
irf_a12 <- irf(var_a4, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # multivariate model

# pre-lisbon period

## response: time pressure invocations; impulse: asylum numbers
irf_a13 <- irf(var_a5, impulse="ts_asylum", response="ts_response", ortho = FALSE,
               cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_a14 <- irf(var_a8, impulse="ts_asylum", response="ts_response", ortho = FALSE,
               cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: time pressure invocations; impulse: mip
irf_a15 <- irf(var_a6, impulse="ts_mip", response="ts_response", ortho = FALSE,
               cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_a16 <- irf(var_a8, impulse="ts_mip", response="ts_response", ortho = FALSE,
               cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: asylum numbers; impulse: time pressure invocations
irf_a17 <- irf(var_a5, impulse="ts_response", response="ts_asylum", ortho = FALSE,
               cumulative = TRUE, seed = 1) # bivariate model
irf_a18 <- irf(var_a8, impulse="ts_response", response="ts_asylum", ortho = FALSE,
               cumulative = TRUE, seed = 1) # multivariate model

## response: asylum numbers; impulse: mip
irf_a19 <- irf(var_a8, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
               cumulative = TRUE, seed = 1) # bivariate model
irf_a20 <- irf(var_a8, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
               cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: time pressure invocations
irf_a21 <- irf(var_a6, impulse="ts_response", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # bivariate model
irf_a22 <- irf(var_a8, impulse="ts_response", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: asylum numbers
irf_a23 <- irf(var_a7, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # bivariate model
irf_a24 <- irf(var_a8, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # multivariate model

# post-lisbon period

## response: time pressure invocations; impulse: asylum numbers
irf_a25 <- irf(var_a9, impulse="ts_asylum", response="ts_response", ortho = FALSE,
               cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_a26 <- irf(var_a12, impulse="ts_asylum", response="ts_response", ortho = FALSE,
               cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: time pressure invocations; impulse: mip
irf_a27 <- irf(var_a10, impulse="ts_mip", response="ts_response", ortho = FALSE,
               cumulative = TRUE, ci = 0.95, seed = 1) # bivariate model
irf_a28 <- irf(var_a12, impulse="ts_mip", response="ts_response", ortho = FALSE,
               cumulative = TRUE, ci = 0.95, seed = 1) # multivariate model

## response: asylum numbers; impulse: time pressure invocations
irf_a29 <- irf(var_a9, impulse="ts_response", response="ts_asylum", ortho = FALSE,
               cumulative = TRUE, seed = 1) # bivariate model
irf_a30 <- irf(var_a12, impulse="ts_response", response="ts_asylum", ortho = FALSE,
               cumulative = TRUE, seed = 1) # multivariate model

## response: asylum numbers; impulse: mip
irf_a31 <- irf(var_a11, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
               cumulative = TRUE, seed = 1) # bivariate model
irf_a32 <- irf(var_a12, impulse="ts_mip", response="ts_asylum", ortho = FALSE,
               cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: time pressure invocations
irf_a33 <- irf(var_a10, impulse="ts_response", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # bivariate model
irf_a34 <- irf(var_a12, impulse="ts_response", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # multivariate model

## response: mip; impulse: asylum numbers
irf_a35 <- irf(var_a11, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # bivariate model
irf_a36 <- irf(var_a12, impulse="ts_asylum", response="ts_mip", ortho = FALSE,
               cumulative = TRUE, seed = 1) # multivariate model

```

\newpage

```{r, vec - irf graphics, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# data frame
irf_df_vec <- data.frame(
  vec_model = c("bivariate", "multivariate", "bivariate", "multivariate", 
                "bivariate", "multivariate", "bivariate", "multivariate",
                "bivariate", "multivariate", "bivariate", "multivariate",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)",
                "bivariate (post-lisbon)", "multivariate (post-lisbon)"),
  response = c(irf_a1$response, irf_a2$response, irf_a3$response, irf_a4$response,
               irf_a5$response, irf_a6$response, irf_a7$response, irf_a8$response,
               irf_a9$response, irf_a10$response, irf_a11$response, irf_a12$response,
               irf_a13$response, irf_a14$response, irf_a15$response, irf_a16$response,
               irf_a17$response, irf_a18$response, irf_a19$response, irf_a20$response,
               irf_a21$response, irf_a22$response, irf_a23$response, irf_a24$response,
               irf_a25$response, irf_a26$response, irf_a27$response, irf_a28$response,
               irf_a29$response, irf_a30$response, irf_a31$response, irf_a32$response,
               irf_a33$response, irf_a34$response, irf_a35$response, irf_a36$response),
  impulse = c(irf_a1$impulse, irf_a2$impulse, irf_a3$impulse, irf_a4$impulse,
              irf_a5$impulse, irf_a6$impulse, irf_a7$impulse, irf_a8$impulse,
              irf_a9$impulse, irf_a10$impulse, irf_a11$impulse, irf_a12$impulse,
              irf_a13$impulse, irf_a14$impulse, irf_a15$impulse, irf_a16$impulse,
              irf_a17$impulse, irf_a18$impulse, irf_a19$impulse, irf_a20$impulse,
              irf_a21$impulse, irf_a22$impulse, irf_a23$impulse, irf_a24$impulse,
              irf_a25$impulse, irf_a26$impulse, irf_a27$impulse, irf_a28$impulse,
              irf_a29$impulse, irf_a30$impulse, irf_a31$impulse, irf_a32$impulse,
              irf_a33$impulse, irf_a34$impulse, irf_a35$impulse, irf_a36$impulse),
  irf = c(tail(unlist(irf_a1$irf), n = 1), tail(unlist(irf_a2$irf), n = 1),
          tail(unlist(irf_a3$irf), n = 1), tail(unlist(irf_a4$irf), n = 1),
          tail(unlist(irf_a5$irf), n = 1), tail(unlist(irf_a6$irf), n = 1),
          tail(unlist(irf_a7$irf), n = 1), tail(unlist(irf_a8$irf), n = 1),
          tail(unlist(irf_a9$irf), n = 1), tail(unlist(irf_a10$irf), n = 1),
          tail(unlist(irf_a11$irf), n = 1), tail(unlist(irf_a12$irf), n = 1),
          tail(unlist(irf_a13$irf), n = 1), tail(unlist(irf_a14$irf), n = 1),
          tail(unlist(irf_a15$irf), n = 1), tail(unlist(irf_a16$irf), n = 1),
          tail(unlist(irf_a17$irf), n = 1), tail(unlist(irf_a18$irf), n = 1),
          tail(unlist(irf_a19$irf), n = 1), tail(unlist(irf_a20$irf), n = 1),
          tail(unlist(irf_a21$irf), n = 1), tail(unlist(irf_a22$irf), n = 1),
          tail(unlist(irf_a23$irf), n = 1), tail(unlist(irf_a24$irf), n = 1),
          tail(unlist(irf_a25$irf), n = 1), tail(unlist(irf_a26$irf), n = 1),
          tail(unlist(irf_a27$irf), n = 1), tail(unlist(irf_a28$irf), n = 1),
          tail(unlist(irf_a29$irf), n = 1), tail(unlist(irf_a30$irf), n = 1),
          tail(unlist(irf_a31$irf), n = 1), tail(unlist(irf_a32$irf), n = 1),
          tail(unlist(irf_a33$irf), n = 1), tail(unlist(irf_a34$irf), n = 1),
          tail(unlist(irf_a35$irf), n = 1), tail(unlist(irf_a36$irf), n = 1)),
  lower = c(tail(unlist(irf_a1$Lower), n = 1), tail(unlist(irf_a2$Lower), n = 1),
            tail(unlist(irf_a3$Lower), n = 1), tail(unlist(irf_a4$Lower), n = 1),
            tail(unlist(irf_a5$Lower), n = 1), tail(unlist(irf_a6$Lower), n = 1),
            tail(unlist(irf_a7$Lower), n = 1), tail(unlist(irf_a8$Lower), n = 1),
            tail(unlist(irf_a9$Lower), n = 1), tail(unlist(irf_a10$Lower), n = 1),
            tail(unlist(irf_a11$Lower), n = 1), tail(unlist(irf_a12$Lower), n = 1),
            tail(unlist(irf_a13$Lower), n = 1), tail(unlist(irf_a14$Lower), n = 1),
            tail(unlist(irf_a15$Lower), n = 1), tail(unlist(irf_a16$Lower), n = 1),
            tail(unlist(irf_a17$Lower), n = 1), tail(unlist(irf_a18$Lower), n = 1),
            tail(unlist(irf_a19$Lower), n = 1), tail(unlist(irf_a20$Lower), n = 1),
            tail(unlist(irf_a21$Lower), n = 1), tail(unlist(irf_a22$Lower), n = 1),
            tail(unlist(irf_a23$Lower), n = 1), tail(unlist(irf_a24$Lower), n = 1),
            tail(unlist(irf_a25$Lower), n = 1), tail(unlist(irf_a26$Lower), n = 1),
            tail(unlist(irf_a27$Lower), n = 1), tail(unlist(irf_a28$Lower), n = 1),
            tail(unlist(irf_a29$Lower), n = 1), tail(unlist(irf_a30$Lower), n = 1),
            tail(unlist(irf_a31$Lower), n = 1), tail(unlist(irf_a32$Lower), n = 1),
            tail(unlist(irf_a33$Lower), n = 1), tail(unlist(irf_a34$Lower), n = 1),
            tail(unlist(irf_a35$Lower), n = 1), tail(unlist(irf_a36$Lower), n = 1)),
  upper = c(tail(unlist(irf_a1$Upper), n = 1), tail(unlist(irf_a2$Upper), n = 1),
            tail(unlist(irf_a3$Upper), n = 1), tail(unlist(irf_a4$Upper), n = 1),
            tail(unlist(irf_a5$Upper), n = 1), tail(unlist(irf_a6$Upper), n = 1),
            tail(unlist(irf_a7$Upper), n = 1), tail(unlist(irf_a8$Upper), n = 1),
            tail(unlist(irf_a9$Upper), n = 1), tail(unlist(irf_a10$Upper), n = 1),
            tail(unlist(irf_a11$Upper), n = 1), tail(unlist(irf_a12$Upper), n = 1),
            tail(unlist(irf_a13$Upper), n = 1), tail(unlist(irf_a14$Upper), n = 1),
            tail(unlist(irf_a15$Upper), n = 1), tail(unlist(irf_a16$Upper), n = 1),
            tail(unlist(irf_a17$Upper), n = 1), tail(unlist(irf_a18$Upper), n = 1),
            tail(unlist(irf_a19$Upper), n = 1), tail(unlist(irf_a20$Upper), n = 1),
            tail(unlist(irf_a21$Upper), n = 1), tail(unlist(irf_a22$Upper), n = 1),
            tail(unlist(irf_a23$Upper), n = 1), tail(unlist(irf_a24$Upper), n = 1),
            tail(unlist(irf_a25$Upper), n = 1), tail(unlist(irf_a26$Upper), n = 1),
            tail(unlist(irf_a27$Upper), n = 1), tail(unlist(irf_a28$Upper), n = 1),
            tail(unlist(irf_a29$Upper), n = 1), tail(unlist(irf_a30$Upper), n = 1),
            tail(unlist(irf_a31$Upper), n = 1), tail(unlist(irf_a32$Upper), n = 1),
            tail(unlist(irf_a33$Upper), n = 1), tail(unlist(irf_a34$Upper), n = 1),
            tail(unlist(irf_a35$Upper), n = 1), tail(unlist(irf_a36$Upper), n = 1))
  ) %>%
  mutate(response = recode(response, 
                           "ts_response" = "Time-pressure invocations",
                           "ts_asylum" = "Asylum numbers",
                           "ts_mip" = "Public pressure")) %>%
  mutate(impulse = recode(impulse, 
                          "ts_response" = "Time-pressure invocations",
                          "ts_asylum" = "Asylum numbers",
                          "ts_mip" = "Public pressure")) %>%
  mutate(vec_model =  fct_relevel(vec_model,
                                  "bivariate", "multivariate",
                                  "bivariate (pre-lisbon)", "multivariate (pre-lisbon)",
                                  "bivariate (post-lisbon)", "multivariate (post-lisbon)"))

# plots
ggplot(irf_df_vec %>% filter(response == "Time-pressure invocations"), 
       aes(x = fct_rev(impulse), y = irf, ymin = lower, ymax = upper, colour = vec_model)) +
  geom_linerange(size = 2.5, position = position_dodge(width = -0.55), alpha = 0.4) +  
  geom_point(position = position_dodge(width = -0.55), size = 3.5) +
  geom_hline(aes(yintercept = 0), color = "black", linetype = 2, alpha = 0.75) +
  scale_color_manual(values=c("#000000", "#999999", "#000099", "#0099FF", "#FFCC00", "#FFFF33")) +
  coord_flip() +
  #scale_y_continuous(limits = c(-2, 6), breaks = c(-2, 0, 2, 4, 6)) +
  labs (subtitle = "",
        title = "",
        y = "Predicted time pressure invocations (12 months cumulative)", x = "",
        colour = "VEC Model") +
  theme_bw()

ggplot(irf_df_vec %>% filter(response == "Asylum numbers"), 
       aes(x = fct_rev(impulse), y = irf, ymin = lower, ymax = upper, colour = vec_model)) +
  geom_linerange(size = 2.5, position = position_dodge(width = -0.55), alpha = 0.4) +  
  geom_point(position = position_dodge(width = -0.55), size = 3.5) +
  geom_hline(aes(yintercept = 0), color = "black", linetype = 2, alpha = 0.75) +
  scale_x_discrete(
    labels = c("Time pressure\ninvocations", "Public pressure")
  ) +
  scale_color_manual(values=c("#000000", "#999999", "#000099", "#0099FF", "#FFCC00", "#FFFF33")) +
  coord_flip() +
  #scale_y_continuous(limits = c(-2, 6), breaks = c(-2, 0, 2, 4, 6)) +
  labs (subtitle = "",
        title = "",
        y = "Predicted asylum numbers (12 months cumulative)", x = "",
        colour = "VEC Model") +
  theme_bw()

ggplot(irf_df_vec %>% filter(response == "Public pressure"), 
       aes(x = fct_rev(impulse), y = irf, ymin = lower, ymax = upper, colour = vec_model)) +
  geom_linerange(size = 2.5, position = position_dodge(width = -0.55), alpha = 0.4) +  
  geom_point(position = position_dodge(width = -0.55), size = 3.5) +
  geom_hline(aes(yintercept = 0), color = "black", linetype = 2, alpha = 0.75) +
  scale_x_discrete(
    labels = c("Time pressure\ninvocations", "Asylum numbers")
  ) +
  scale_color_manual(values=c("#000000", "#999999", "#000099", "#0099FF", "#FFCC00", "#FFFF33")) +
  coord_flip() +
  #scale_y_continuous(limits = c(-2, 6), breaks = c(-2, 0, 2, 4, 6)) +
  labs (subtitle = "",
        title = "",
        y = "Predicted public pressure (12 months cumulative)", x = "",
        colour = "VEC Model") +
  theme_bw()

```

\newpage

## A.3 Alternative regression model specifications

The following analysis reports OLS models with different lags to provide additional robustness checks for the findings from the VAR models (see main article) and VECM models (see Appendix A.2). The models investigate the (lagged) effect of asylum numbers and public pressure on the invocation of time pressure from the European Commission.

The first set of models is based on simple OLS regression without lagged variables, while the following sets of models include 1-month, 2-month, 3-month, 6-month and 12-month lags for the main independent variables (i.e. asylum numbers and public pressure) as well as corresponding lagged dependent variables. These specifications allow to test the robustness of the findings across different time lag structures.

Overall, the findings are largely stable across all the specifications and align with the findings from the VAR and VECM models. We find a robust effect of asylum numbers and public pressure on time pressure invocations (see first and second models in Tables A?-A?). These findings also largely hold in models that include both independent variables (see third models in Tables A?-A?, but see fourth model in Table A? and fifth model in Table A?). Furthermore, the models confirm the moderating factor of legal authority on policy-making in the EU. We find a consistent positive interaction effect between the post-Lisbon period and asylum numbers and public pressure for time pressure invocations (see fourth and fifth models in Tables A?-A?).

### OLS

```{r, ols - models, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# ols

## models
ols_1 <- lm(total_response_count ~ 
              asylum_total_scaled + 
              post_lisbon,
            data = main_data)
ols_2 <- lm(total_response_count ~ 
              avg_mip_country_immigration +
              post_lisbon,
            data = main_data)
ols_3 <- lm(total_response_count ~ 
              asylum_total_scaled +
              avg_mip_country_immigration +
              post_lisbon,
            data = main_data)
ols_4 <- lm(total_response_count ~ 
              asylum_total_scaled * post_lisbon,
            data = main_data)
ols_5 <- lm(total_response_count ~ 
              avg_mip_country_immigration * post_lisbon,
            data = main_data)

## output
modelsummary(models = list(ols_1, ols_2, ols_3, ols_4, ols_5),
             coef_map = c("asylum_total_scaled" = "Asylum numbers", 
                          "avg_mip_country_immigration" = "Public pressure",
                          "post_lisbon" = "Post-Lisbon",
                          "asylum_total_scaled:post_lisbon" = 
                            "Asylum numbers:Post-Lisbon",
                          "avg_mip_country_immigration:post_lisbon" = 
                            "Public pressure:Post-Lisbon",
                          "(Intercept)" = "(Intercept)"),
             output = "kableExtra", 
             stars = TRUE) %>%
  add_header_above(c(" " = 1, "DV: Time pressure invocations" = 5)) %>%
  kable_styling(latex_options = c("scale_down", "HOLD_position"))

```

\newpage

### OLS (lag -1)

```{r, ols lag1 - models, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# ols (lag -1)

## models
ols_lags_1 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag1 +
                   post_lisbon +
                   total_response_count_lag1,
                 data = main_data)
ols_lags_2 <- lm(total_response_count ~ 
                   avg_mip_country_immigration_lag1 + 
                   post_lisbon +
                   total_response_count_lag1,
                 data = main_data)
ols_lags_3 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag1 + 
                   avg_mip_country_immigration_lag1 +
                   post_lisbon +
                   total_response_count_lag1,
                 data = main_data)
ols_lags_4 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag1 * post_lisbon + 
                   total_response_count_lag1,
                 data = main_data)
ols_lags_5 <- lm(total_response_count ~ 
                   avg_mip_country_immigration_lag1 * post_lisbon +
                   total_response_count_lag1,
                 data = main_data)

## output
modelsummary(models = list(ols_lags_1, ols_lags_2, ols_lags_3, ols_lags_4, ols_lags_5),
             coef_map = c("asylum_total_scaled_lag1" = "Asylum numbers (t-1)", 
                          "avg_mip_country_immigration_lag1" = "Public pressure (t-1)",
                          "post_lisbon" = "Post-Lisbon",
                          "asylum_total_scaled_lag1:post_lisbon" = 
                            "Asylum numbers (t-1):Post-Lisbon",
                          "avg_mip_country_immigration_lag1:post_lisbon" = 
                            "Public pressure (t-1):Post-Lisbon",
                          "total_response_count_lag1" = 
                            "Time pressure invocations (t-1)",
                          "(Intercept)" = "(Intercept)"),
             output = "kableExtra", 
             stars = TRUE) %>%
  add_header_above(c(" " = 1, "DV: Time pressure invocations" = 5)) %>%
  kable_styling(latex_options = c("scale_down", "HOLD_position"))

```

\newpage

### OLS (lag -2)

```{r, ols lag2 - models, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# ols (lag -2)

## models
ols_lags_1 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag2 +
                   post_lisbon +
                   total_response_count_lag2,
                 data = main_data)
ols_lags_2 <- lm(total_response_count ~ 
                   avg_mip_country_immigration_lag2 + 
                   post_lisbon +
                   total_response_count_lag2,
                 data = main_data)
ols_lags_3 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag2 + 
                   avg_mip_country_immigration_lag2 +
                   post_lisbon +
                   total_response_count_lag2,
                 data = main_data)
ols_lags_4 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag2 * post_lisbon + 
                   total_response_count_lag2,
                 data = main_data)
ols_lags_5 <- lm(total_response_count ~ 
                   avg_mip_country_immigration_lag2 * post_lisbon +
                   total_response_count_lag2,
                 data = main_data)

## output
modelsummary(models = list(ols_lags_1, ols_lags_2, ols_lags_3, ols_lags_4, ols_lags_5),
             coef_map = c("asylum_total_scaled_lag2" = "Asylum numbers (t-2)", 
                          "avg_mip_country_immigration_lag2" = "Public pressure (t-2)",
                          "post_lisbon" = "Post-Lisbon",
                          "asylum_total_scaled_lag2:post_lisbon" = 
                            "Asylum numbers (t-2):Post-Lisbon",
                          "avg_mip_country_immigration_lag2:post_lisbon" =
                            "Public pressure (t-2):Post-Lisbon",
                          "total_response_count_lag2" = 
                            "Time pressure invocations (t-2)",
                          "(Intercept)" = "(Intercept)"),
             output = "kableExtra", 
             stars = TRUE) %>%
  add_header_above(c(" " = 1, "DV: Time pressure invocations" = 5)) %>%
  kable_styling(latex_options = c("scale_down", "HOLD_position"))

```

\newpage

### OLS (lag -3)

```{r, ols lag3 - models, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# ols (lag -3)

## models
ols_lags_1 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag3 +
                   post_lisbon +
                   total_response_count_lag3,
                 data = main_data)
ols_lags_2 <- lm(total_response_count ~ 
                   avg_mip_country_immigration_lag3 + 
                   post_lisbon +
                   total_response_count_lag3,
                 data = main_data)
ols_lags_3 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag3 + 
                   avg_mip_country_immigration_lag3 +
                   post_lisbon +
                   total_response_count_lag3,
                 data = main_data)
ols_lags_4 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag3 * post_lisbon + 
                   total_response_count_lag3,
                 data = main_data)
ols_lags_5 <- lm(total_response_count ~ 
                   avg_mip_country_immigration_lag3 * post_lisbon +
                   total_response_count_lag3,
                 data = main_data)

## output
modelsummary(models = list(ols_lags_1, ols_lags_2, ols_lags_3, ols_lags_4, ols_lags_5),
             coef_map = c("asylum_total_scaled_lag3" = "Asylum numbers (t-3)", 
                          "avg_mip_country_immigration_lag3" = "Public pressure (t-3)",
                          "post_lisbon" = "Post-Lisbon",
                          "asylum_total_scaled_lag3:post_lisbon" = 
                            "Asylum numbers (t-3):Post-Lisbon",
                          "avg_mip_country_immigration_lag3:post_lisbon" = 
                            "Public pressure (t-3):Post-Lisbon",
                          "total_response_count_lag3" = 
                            "Time pressure invocations (t-3)",
                          "(Intercept)" = "(Intercept)"),
             output = "kableExtra", 
             stars = TRUE) %>%
  add_header_above(c(" " = 1, "DV: Time pressure invocations" = 5)) %>%
  kable_styling(latex_options = c("scale_down", "HOLD_position"))

```

\newpage

### OLS (lag -6)

```{r, ols lag6 - models, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# ols (lag -6)

## models
ols_lags_1 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag6 +
                   post_lisbon +
                   total_response_count_lag6,
                 data = main_data)
ols_lags_2 <- lm(total_response_count ~ 
                   avg_mip_country_immigration_lag6 + 
                   post_lisbon +
                   total_response_count_lag6,
                 data = main_data)
ols_lags_3 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag6 + 
                   avg_mip_country_immigration_lag6 +
                   post_lisbon +
                   total_response_count_lag6,
                 data = main_data)
ols_lags_4 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag6 * post_lisbon + 
                   total_response_count_lag6,
                 data = main_data)
ols_lags_5 <- lm(total_response_count ~ 
                   avg_mip_country_immigration_lag6 * post_lisbon +
                   total_response_count_lag6,
                 data = main_data)

## output
modelsummary(models = list(ols_lags_1, ols_lags_2, ols_lags_3, ols_lags_4, ols_lags_5),
             coef_map = c("asylum_total_scaled_lag6" = "Asylum numbers (t-6)", 
                          "avg_mip_country_immigration_lag6" = "Public pressure (t-6)",
                          "post_lisbon" = "Post-Lisbon",
                          "asylum_total_scaled_lag6:post_lisbon" =
                            "Asylum numbers (t-6):Post-Lisbon",
                          "avg_mip_country_immigration_lag6:post_lisbon" = 
                            "Public pressure (t-6):Post-Lisbon",
                          "total_response_count_lag6" = 
                            "Time pressure invocations (t-6)",
                          "(Intercept)" = "(Intercept)"),
             output = "kableExtra", 
             stars = TRUE) %>%
  add_header_above(c(" " = 1, "DV: Time pressure invocations" = 5)) %>%
  kable_styling(latex_options = c("scale_down", "HOLD_position"))

```

\newpage

### OLS (lag -12)

```{r, ols lag12 - models, include=TRUE, echo=TRUE, warning=FALSE, message=FALSE, dev=c('pdf', 'png')}

# ols (lag -12)

## models
ols_lags_1 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag12 +
                   post_lisbon +
                   total_response_count_lag12,
                 data = main_data)
ols_lags_2 <- lm(total_response_count ~ 
                   avg_mip_country_immigration_lag12 + 
                   post_lisbon +
                   total_response_count_lag12,
                 data = main_data)
ols_lags_3 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag12 + 
                   avg_mip_country_immigration_lag12 +
                   post_lisbon +
                   total_response_count_lag12,
                 data = main_data)
ols_lags_4 <- lm(total_response_count ~ 
                   asylum_total_scaled_lag12 * post_lisbon + 
                   total_response_count_lag12,
                 data = main_data)
ols_lags_5 <- lm(total_response_count ~ 
                   avg_mip_country_immigration_lag12 * post_lisbon +
                   total_response_count_lag12,
                 data = main_data)

## output
modelsummary(models = list(ols_lags_1, ols_lags_2, ols_lags_3, ols_lags_4, ols_lags_5),
             coef_map = c("asylum_total_scaled_lag12" = "Asylum numbers (t-12)", 
                          "avg_mip_country_immigration_lag12" = "Public pressure (t-12)",
                          "post_lisbon" = "Post-Lisbon",
                          "asylum_total_scaled_lag12:post_lisbon" = 
                            "Asylum numbers (t-12):Post-Lisbon",
                          "avg_mip_country_immigration_lag12:post_lisbon" = 
                            "Public pressure (t-12):Post-Lisbon",
                          "total_response_count_lag12" = 
                            "Time pressure invocations (t-12)",
                          "(Intercept)" = "(Intercept)"),
             output = "kableExtra", 
             stars = TRUE) %>%
  add_header_above(c(" " = 1, "DV: Time pressure invocations" = 5)) %>%
  kable_styling(latex_options = c("scale_down", "HOLD_position"))

```
